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j* 1 Abstract 
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This work addresses linear transport in turbulent media, with emphasis on neutral particle 

e 

(atoms, molecules) transport in magnetized fusion plasmas. A stochastic model for turbulent 
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plasmas, based upon a multivariate Gamma distribution, is presented. The geometry is a 2D 

slab and turbulence is assumed to be statistically homogeneous. The average neutral density 
• i— i 

and ionization source, which are the quantities relevant for integrated simulations and diagnostic 
applications, are calculated analytically in the scattering free case. The boundary conditions and 
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the ratio of the turbulence correlation length to the neutral mean free path are identified as the 
main control parameters in the problem. The non trivial relationship between the average neutral 
density and the ionization source is investigated. Monte Carlo calculations including scattering are 
then presented, and the main trends obtained in the scattering free case are shown to be conserved. 
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I. INTRODUCTION 

The problem of linear transport in a stochastic background medium has attracted much 
interest in the last twenty years, with applications ranging from neutronics in boiling fluids 
to radiation transport in cloudy atmospheres PQ . Technically, the problem consists in treat- 
ing the absorption and scattering rates, which describe the interaction of particles with the 
host medium in the linear transport equation, as random variables. The resulting stochastic 
integro-differential equation is then solved for the statistical moments of the particle distri- 
bution function. The back-reaction of transport on the host medium is neglected, and the 
statistical properties of the latter have to be specified for each particular problem. The most 
thoroughly studied case is a mixture of two fluids, where the rates are discrete two states 
random variables pp. The case of Gaussian statistics has subsequently been studied, as a 
first attempt to address the problem of neutral particle transport in fusion plasmas [2j [3] , 
where plasma turbulence leads to a continuous probability density function for the rates. In 
fact, in magnetic fusion devices such as tokamaks, the transport of species such as atoms 
or molecules, as well as line radiation, is best described in kinetic terms. The single par- 
ticle distribution function /(v,r,t) obeys a Boltzmann equation, which is linear provided 
neutral-neutral collisions are negligible. In practice, given the complexity of the geome- 
try and the numerous reaction channels involved, Boltzmann's equation is often solved by a 
Monte Carlo approach [I] . Several computer codes have been developed also in the magnetic 
fusion context, among which the EIRENE code [5] which is used for ITER modeling, often 
iteratively in conjunction with the B2 plasma fluid code [6]. Indeed, at the outer plasma 
region (the so called "plasma edge" ) , where the density of neutrals can be of the same order 
as the electron density, the fluid equations describing the plasma and the kinetic neutral 
transport problem must be solved consistently. These models have a mesoscopic character, 
in the sense that the time scales of interest (of the order of 10 ms) are large compared 
to time scales characterizing turbulence (correlation time of about 10 fis [7]), and small 
compared to macroscopic times in the discharge (e.g. seconds). The effect of turbulence at 
these mesoscopic scales is essential, since turbulent transport usually dominates collisional 
transport [8j. Turbulent transport is described in current fluid codes by ad-hoc empirical 
transport coefficients (Dj_, v±) in the direction transverse to the magnetic field, such that for 
the particle flux T turb = (Nv) = —D±V(N) +v±(N), where brackets denote a time average 



(N and v stand for the fluctuating part of the density and the velocity fields). However, 
neutral transport is still commonly calculated on the average plasma background, without 
any account of the underlying turbulence. The problem is especially acute at the plasma 
edge, since edge plasma turbulence has different properties compared to core turbulence 
(prevailing in the central part of magnetically confined plasmas). In particular, fluctuation 
rates up to order unity are observed in the outer edge region of the plasma, where magnetic 
field lines are "open" (scrape off layer, SOL), i.e. intersect solid material components. These 
large fluctuation rates are often interpreted in terms of propagating plasma filaments (aka 
"blobs" or "avalanches"), which have been seen both experimentally [9l [TTJ] and numerically 
|llj . In the drift plane, perpendicular to the magnetic field, these filaments have an ex- 
tension of about 1 cm. They propagate radially outward with a velocity of the order of 1 
km/s [7]. The question of the existence of a significant contribution of turbulence to neutral 
particle transport thus naturally arises. Moreover, plasma spectroscopy usually works with 
integration times much larger than the turbulent time scales, so that the measured signals, 
such as Doppler broadened spectral line profiles, are time and space averages over fluctua- 
tions [12]. As a result, the observed Doppler profile is related to the average kinetic velocity 
distribution (/). In the following, the time averages involved in either transport or spectro- 
scopic data modeling will be replaced by an ensemble average over a number of realizations 
of the turbulent density fields. The most attractive feature of such a stochastic model is 
its flexibility, which allows to fully explore the physics of the problem, and in particular to 
study limiting cases. However, the model parameters can also be specified from experimen- 
tal and/or numerical results. The outline of the work is the following. Section [H] is devoted 



to setting up the problem of kinetic transport in a stochastic background. In section HI 
the stochastic model used to describe the turbulent background medium is presented. This 
model, which relies on a multivariate Gamma distribution, is consistent with experimental 
data and alleviates difficulties encountered with Gaussian statistics. Section [IV] deals with 
the numerical implementation of the model in the EIRENE Monte Carlo solver, and its 
validation against analytical results. Finally numerical calculations including scattering are 
presented. 



II. STOCHASTIC LINEAR BOLTZMANN EQUATION 

A. General formulation 

The transport of neutral particles in a fusion plasma is best described in terms of a 
Boltzmann equation 

| + ^.v)Ar,«,n,«)=(|) e , (i) 

where /(r, v, f2, t) is the neutral particles velocity distribution at point r, v the velocity and 
f2 the unit vector along v, such that v = vfl. Inhomogeneous terms due to primary sources 
of neutral particles, are omitted here and in the following discussions, for simplicity. For 
ionizing conditions, the collision operator on the r.h.s. of the previous equation is given for 
hydrogen isotopes atoms by 



§0 = "Oio + v<*)f + f°° v' 2 dv' J dnV CT (|v - v'|)|v - v'|/(r, v', tf)fi(r, v, n), (2) 

where z/j and v cx are respectively the ionization (by electron impact) and charge-exchange 
rates, a cx is the charge exchange cross section, and fi the ion velocity distribution. The 



rates are defined by v io = N e <Ji v e and v cx = Ni(j cx \v — v%\ where the overbars respectively 
stand for an average over the electron and the ion velocity distributions (assumed to be 
maxwellian), and Oi is the ionization cross section. For H 2 (and D 2 , DH, ...) molecules, a 
reasonable first approximation is 
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where Ud is an effective dissociation rate, including contributions from ionization to H 2 ■, dis- 
sociation by electron impact, dissociative ionization, possibly ion conversion (see [13] for an 
overview and more advanced models). A source term describing atoms formed by molecular 
dissociative processes should be added to Eq. rt2J) to treat both species simultaneously, but 
this will not be addressed analytically. Eq. ([T| together with Eq. ([2]) (or Eq. ^) form 
the basis of the linear transport problem [Til [15] (linear in the sense that the collision term 
given by Eq. (J2J) and (J3J) are linear in /). In a general theoretical setting, ionization and 
charge exchange would respectively be called absorption and scattering. Although in the 



following we discuss neutral particle transport in turbulent plasmas, most of the results are 
presented in normalized units, so as to be independent on the precise value of the rates. 
Therefore, these results could be of general interest for linear transport theory. 

The rate coefficients and rates (frequencies) are known functions of the plasma density 
and temperature, and in case of charge exchange also of the neutral velocity v. When 
the plasma is turbulent, the values of the rates also reflect the plasma fluctuations, thus 
introducing a coupling between plasma turbulence and neutral particle transport. In this 
work, we will neglect the back reaction of the resulting fluctuations of the neutral density or 
temperature on the plasma turbulence itself, i. e. neutral particles will be treated as a passive 
species. The statistical properties of turbulence are then input quantities, determined by 
independent turbulence modeling or inferred from experiments. In the following, we will 
assume that the problem is stationary in time, that is the time derivative will be neglected 
in Eq. (II]). This approximation is justified whenever the typical life time of a neutral 
particle, given by z/" 1 , is much smaller than the typical time over which turbulent fields 
evolve (a few microseconds). This is in general not the case in fusion edge plasmas, but this 
approximation is very useful to understand the physics. It could be relaxed in the numerical 



approach presented in section |IV[ which relies on the Monte Carlo code EIRENE. For the 
time being, we are interested in the time average of the neutral distribution over durations 
much larger than the turbulence correlation time, so that according to the ergodic theorem, 
time averages can be replaced by ensemble averages provided the statistics is stationary [16J. 
Averaging Eq. ([!]) over turbulent fluctuations yields 

v ,2 dv' J dSl'a cx (\v-v'\)\v-v'\{f(r,v',Sl , )f i (r,v,Sl)). 

(4) 
This is clearly not a closed equation for (/), since the averages on the r.h.s are not given as 
functions of (/). Neglecting turbulence, as is currently common practice in neutral particle 
transport modeling for tokamaks {e.g. in all B2-EIRENE, EDGE-2D-EIRENE or EMC3- 
EIRENE applications) for neutrals, amounts to solve Boltzmann's equation using the value 
of the rates calculated for the average values of the turbulent fields. For instance, the first 
term on the r.h.s. of Eq. Q is in fact currently approximated by 



{iyio + v cx )f) ^ K((N e ), (T e )) + u cx ((N t ), (Ti))] (f). (5) 

Assuming stationarity, Eq. (fil) is equivalent to the following integral equation 
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where /o = f(0,v,fl), and Q cx is the source of neutrals with velocity v in direction Q at 
point r due to charge exchange, given by 



Q cx (r,v,fl) 
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fi(r,v,Q), (7) 



and s the distance to the boundary along direction i7, i.e. s = |r — rj,| (where r& is the 
intersection between the neutral particle trajectory and the boundary, see Fig. [l|. t(s) is 
the optical thickness defined by 
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where i/ = Ui + z/ cx is the total attenuation rate. The exponential terms in Eq. (pi) give 
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FIG. 1: Sketch of the different vectors involved in the integral formulation of the problem, r defines 
the position at which the neutral particle distribution is evaluated, rb the position of the boundary 
along the neutral flight assuming their velocity is parallel to the unit vector fl. s' measures the 
distance of the neutral flight, starting from r, i.e. going backwards along the neutral particle 
trajectory. 



the probability that a neutral particle does not undergo a collision during a flight at least 



of length s, so that Eq. (JsT) has a natural probabilistic interpretation. Finally, Eq. Q is 
equivalent to the following integral equation 

P s rln' i \ 

(f(r,v,n)) = (f e^))+ / —/ Qcx{r -s>n,v,n)e-^), (9) 
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obtained by averaging Eq. ^ over the background fluctuations. The average neutral density 

is defined by 
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(JV(r))= / v 2 dv / dQ (f(r,v,Q)). (10) 
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In the case of atoms in a plasma, (N(r)) may be determined by Laser Induced Fluorescence 
(LIF) [T7j, by averaging successive measurements performed in a plasma in statistical steady 
state. The primary quantity of interest for integrated edge plasma engineering codes is the 
average ionization source (S(r)) = (u(r)N(r)), which provides the local particle source in 
the plasma equations. Its volume integral is equal to the average of the source flux, thus 
ensuring that the total neutral content (N tot ) is stationary. It should be emphasized that 
(S) is in general NOT related to the average density in a simple and obvious way, in other 
words: (S) ^ (y){N). Instead, (N to t) depends on the properties of the fluctuations. In the 
context of radiation transport, (S) would for instance be an average photo-excitation rate. 

B. Geometry 

We consider a 2D slab geometry, of size C x C in the xy plane. The particle source 
is located at x — 0, and will be specified in section II D| This geometry will be further 



degraded to a ID problem as a first step in the analytical calculations. When dealing with 
neutral particle transport in tokamak plasmas in particular, the x axis represents the radial 
direction, along the minor radius of the torus. The y axis is the poloidal direction, and the z 
direction is along the magnetic field (ignoring the typically small field line pitch). Turbulent 
wave numbers in the parallel direction fen are such that k\\ 3> k±, so that the plasma is taken 
to be homogeneous along z ("flute approximation"), which justifies a 2D model [TTJ. 3D 
effects are believed to be essential to reproduce experimental features of edge turbulence 
[18] . but for our purposes where the statistical properties of turbulence are imposed, the 
2D approximation is sufficient. In the following, our work will focus on the average neutral 
density and ionization source as a function of the radial coordinate x. (N(x)) is related to the 



flux of neutrals crossing a given radial surface x = est. In the simplest case, i.e. one speed 
scattering free transport problem in ID, T(x) = (N(x))vq where vq is the neutral velocity. 
As a result, the screening efficiency p(x) of slab of thickness x, p(x) = 1 — r(x)/T(0), is 
directly linked to (N(x)}. 

C. Stochastic model for turbulent plasma background 

The statistical properties of the rate coefficients are directly related to those of the plasma 
parameters, since for example v = Ui (N e ,T e ). In principle, calculating the averages in Eq. 
(|9| requires knowledge of the functional PDF for the field v(r), W[u]. In practice, this 
functional is not known, but the PDF at one point P\(y, r) {P\{y, r)du gives the probability 
that the rate takes a value between v and v + dv at point r ), and the 2 point correlation 
function C(r, r') = (v(r)u(r')) — (i/(r))(z/(r')) (related to the spatial spectra obtained by a 
Fourier transform) can be measured or calculated numerically In the following, we assume 
statistical homogeneity and isotropy, so that C(r, r') = C(|r — r'|). To avoid unnecessary 
complications related to functional integration, we discretize the problem in space, and focus 
on a ID geometry first. The number of cells n in the spatial grid has to be chosen large 
enough so that the cell size e is much smaller than the turbulence correlation length A. In this 
approximation, W\y\ becomes a multivariate distribution W n (u) where u — (u±, . . . , v n ), v- h 
being the rate value in cell number i. The simplest choice for W n would seem to be a Gaussian 
multivariate distribution, as in references [2J [3] . This has the advantage to permit an exact 
analytical closure of Eq. Q, through an elegant use of Novikov's theorem |19J . The average 
velocity distribution (/) then obeys an effective Bolzmann equation, in which the absorption 
and scattering rates are renormalized by fluctuations. However, the absorption rate becomes 
negative for large fluctuation rates, a consequence of an unphysical feature of the model. In 
fact, for Gaussian statistics, there are realizations for which the rates are negative, at least in 
some regions of space. The weight of these negative values is all the larger that the fluctuation 
rate increases. This problem is very well illustrated in Ref. [20], where it is shown that the 
Gaussian model might nevertheless provide interesting qualitative information on the physics 
of the problem. In summary, the main drawback of the Gaussian model is its limited ability 
to describe the large fluctuation rate regime, where the effects of fluctuations become really 
significant. Truncating the Gaussian distribution to positive variables is not a good solution, 
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because the nice properties of Gaussian statistics are then lost. In particular, even in the 
scattering free case, this choice precludes analytical calculations, since the integrals to be 
calculated are then orthant probabilities ([2T], p. 127). As a result, in the following we 
will consider a statistical model involving only positive random variables. Moreover, recent 
experimental investigations of edge turbulence have shown that density fluctuations in the 
SOL are not Gaussian, and suggest at least two possible choices. In the SOL, the PDF Pi 
has been found to be well fitted by either Log- normal or Gamma distributions (e.g. [22]). 
Taking W n as multivariate Log-normal distribution ^T\ precludes analytical calculations 
of the averages involved in Eq. fcfy. Therefore, for validation purposes, in the following 
we focus on a multivariate generalization of the Gamma distribution (in fact of the chi-2 
distribution), which is the marginal of the Wishart nxn matrix distribution, for which only 
the n diagonal elements are retained [23]. To define the multivariate Gamma distribution, 
first consider the following n-variate Gaussian PDF 

PJX) = \ exp f--X T G" 1 X^) , (11) 

1 ' (27r)"/ 2 ydeT(G) V 2 J 1 K ' 

where G h j = (XiXj) and (JQ) = 0. Let Xy, i — 1, . . . , n and j — 1, . . . , M be M indepen- 
dent series of n zero average Gaussian random numbers sampled from P n (i.e. such that 
(XijXim) = GuSjm), then the n variables 

M 

"* = £*«' (12) 

are distributed according to W n (u). The latter is such that its 1 point marginal W\{yi) 
is a chi-2 distribution with M degrees of freedom, i.e. a Gamma distribution of shape 
parameter /3 = M/2 (see Appendix A), and its 2 point correlation function is given by 



(i v i u j)) = ( u i u j) ~ ( v i)( v j) = 2MGfj. In fact from Eq. (12), we have (ui) = MGu and 
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((^» = Y,i x li x ij) + E<*«><*3> - M2 °^ G n- (13) 

Expressing the fourth moment (X^Xfj) in terms of the second moments, we get {{vk.vi)) = 
2MG 2 kl . In other words, the correlation matrix of the v variables is CV,- = 2MG?-, so that 
it is fully specified by Gij. However, note that G^ is not defined in a unique way from CV,-. 
Moreover, Cy > by construction, so that the multivariate Gamma model cannot describe 



situations where the correlation function oscillates around 0. As an example, to obtain the 
correlation matrix C corresponding to 

C(r-r') =a 2 exp(-|r-r , |/A), (14) 

where a is the standard deviation, and A the turbulence correlation length (or integral scale), 
one should set Cy = C(r$ — r 3 -), where r t and r,,- are the coordinates of the center of the cells 
% and j, respectively. Then, 

^ = ^exp(-|r 2 ,-r;.|/2A). (15) 

The fluctuation rate R (%) is determined by M through R = ■\J2jM x 100, and in the 
following we shall limit ourselves to M > 3, so that W\(vi = 0) = 0. It should be noted that 
an integral scale can be defined for any correlation function by 



A= / C(r)dr/C(0), (16) 
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provided the integral converges. A represents the typical size of the turbulent structures. 
The case where A is infinite, because of long range correlations (e.g. C(r) oc r~ a , a < 1), 
will not be considered here. Instead, the limit A — > +oo should be understood as a situation 
where the plasma background becomes spatially homogeneous, see Fig. [8j The most useful 
form of W n {y) for our purposes is the following, 
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where P n is given by Eq.(ll) . Eq. (17) is just the mathematical translation of the definition 
of the statistics of j/j. The generalization to a 2D domain consisting of n x n cells can be made 
in two equivalent ways, either by using Eq. (17) with n replaced by n 2 and a proper labeling 



scheme for the cells, or by using Wishart's distribution (for which the random variable is a 
n x n matrix) |21j . 

D. Sources 

We first address the question of boundary conditions (i.e. more precisely sources in the 
case where charge exchange is included) from the physics rather than the technical point 

10 



of view, hence focusing on neutral particles in plasma. The main source of neutrals in a 
magnetic fusion confinement device is recycling, that is ions hitting the wall subsequently 
reenter the plasma as neutrals [23]. Recycling is a complex process, involving several different 
elementary mechanisms, the time scales of which are not necessarily short compared to those 
of turbulence. For instance, ions impinging on the wall can be backscattered as neutrals, and 
this fraction of the recycling flux will be directly related to the instantaneous plasma flux 
because time scales involved are very short (< 1CT 9 s). The backscattering probability is of 
the order of 0.1 for hydrogen on carbon material surfaces, and can be significantly larger for 
heavier fusion relevant materials like tungsten [23] ■ Hydrogen ions can, after neutralization, 
also thermalize with the wall, then recombine with another atom in a H2 molecule, which 
then desorbs. In this case, time scales involved depend on the wall properties, i.e. whether it 
is saturated with hydrogen or not. In the latter case, the wall acts as a pump and its response 
becomes very slow [25]. As a result, in the following, we will limit ourselves to two extreme 
cases that will be called slow and fast recycling. In the first case, recycling is slow compared 
to turbulence, so that neutrals are emitted homogeneously from the wall (the average plasma 
flux is assumed to be homogenous along y). In this case, the boundary condition at the wall 
is a non stochastic and uniform flux at the wall. The second case corresponds to a situation 
where recycling is much faster than turbulence, so that the local neutral flux at the wall is 
a direct reflection of the plasma flux distribution T p . The latter is assumed to be given by 
T p (y) = N e (0,y)Vuob (unit m _2 .s _1 ), where Vuob is the typical "plasma blob" velocity, of the 
order of Vuob — 1 km/s [7J. In reality if the wall is not saturated, a superposition of these 
two limiting cases may provide a reasonable description of the recycling source. Finally, it 
should be pointed out that in the frame of radiation transport, the fast recycling case would 
describe a reflective boundary, while the slow recycling case would represent the surface of 
a black body. 

III. ANALYTICAL CASE FOR VALIDATION 

In order to obtain a model which can be solved analytically, we consider the scattering 
free (i.e. no charge exchange) problem in a 2D slab geometry. Only density fluctuations are 
considered, and the ionization rate v is assumed to be linear in density (which physically 
means that the multi-step effects of the collisional radiative ionization cascade are neglected). 
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In this case, the multivariate PDFs of the rate coefficient and of the density are identical 
up to a scaling factor. We obtain the solution in two steps. In the scattering free problem, 
neutral particles move along straight lines, and we start by solving for the neutral velocity 
distribution along a given direction, defined by the unit vector f2 (ID problem). We further 
assume that the source is mono-kinetic, with a velocity v = v o (one speed transport problem). 
In this simplified case, Eq. ^ reduces to 

(f(a,v,a)) = -S(v-Vo)(T e-<^), (18) 



where r is the optical thickness defined by Eq. rt8j). T is the flux density (unit m~ 2 .s -1 ) 
at the wall (x = 0), which is either non-stochastic (slow recycling) or proportional to the 
instantaneous plasma density at the wall (fast recycling). The average neutral particle 



density in 2D (N2d(%)} is obtained by first averaging Eq. (18) on vq and f2, then integrating 
over v. The average ionization source is given by 

(S(x)) = Mx)N(x)) = - (T u(x)e- T(x) ) . (19) 

Vq 

The scattering free case is very interesting both for code validation purposes and for un- 
derstanding the physics of the problem. Moreover, it provides a reasonable approximation 
for neutral species sputtered from the wall, and for H2 molecules for which scattering can 
indeed be neglected, except in high density and low temperature plasmas (N e > 10 20 m -3 
and T e < 5 eV), where elastic collisions have to be retained |13j . 

A. General results 

There are two important results which are valid for any choice of the turbulence statistical 
properties. The first one pertains to the vanishing turbulence correlation length A limit, that 
is when the ratio a of A to the neutral mean free path I tends to zero. In fact, the quantity 
t(s)/s, where r(s) is the optical thickness defined by Eq. (J8|, is nothing else but the ID 
spatial average of u(r — s'Q)/v along the neutral particle trajectory. In the limit where 
X/s — > 0, the ergodic theorem [16] ensures that for almost all realizations of v 

- [ S u((r-s'n)/v )ds' = (r(s)), (20) 

s Jo 
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which implies 

(f(x))x = o = ^S(v-v )e-^l (21) 

v 

Therefore, in the A — >• limit, the average neutral distribution obeys an effective (scattering 
free) Boltzmann equation in which the ionization rate is simply replaced by its average value. 
In other words, the crossed term (i/f) factorizes in (i/)(f). This limit is called the atomic 
mix limit in Ref. jT], where it is obtained by different arguments. From a physical point of 
view, at a distance x ^> A from the source, the correlations between the values of / and v 
are expected to be negligible. Indeed, / depends on the values taken by v between and x. 
It should be noted that in general the ionization rate v is a non linear function of the plasma 
parameters (N e , T e ), so that (u(N e ,T e )) ^ u((N e ), (T e )). Therefore, the turbulence free case 
is recovered in the A — > limit only when v is a linear function of N e (and more generally, 
when temperature fluctuations are neglected). Now, we come to the second important point. 
Jensen's inequality [2S], which pertains to convex functions, states that 

<e" T(x) ) ^ e~ {T{x)) . (22) 



The equality is realized in the A = limit (see Eq. (21)), i.e. the fastest decaying average 
density profile is obtained for vanishing correlation lengths. In the model considered in this 
section (one speed transport, no scattering, v linear in density), we have thus proved that, 
provided temperature fluctuations are neglected, the optical thickness of the plasma t(s) is 
always reduced by fluctuations. 

B. Characteristic function 

The quantities we are interested in, namely the average density (N) and ionization source 
(S) = (vN), are related to the characteristic functional associated to W[v], 

Z[u] = (exp- v{y)u{y)dy\ , (23) 

and to its functional derivatives with respect to u, where u(y) is an arbitrary function of y. 
For instance, 
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8Z[u] 
Su(a) 



(u = Q(x' — x)/v ) 



u(a)exp- / z/(|/) — 
o v o 



(24) 



where B(x' — x) is the Heaviside step function. It turns out that the discretized version of 



Z[u], Z n {y) with v 



\y\, 



being a column vector, can be calculated analytically, 



provided the rate coefficients statistics is given by a multivariate Gamma distribution. Z n 
is defined by 



ZJvl] 



duW n (u) exp v ■ u, 

v 



(25) 



«!,... , u n ) and e is the discretization step. Plugging Eq. (17) and (11) into 



where u 

the previous equation, and inverting the order of integration leads to 



Z n (u) 



1 



J eOCexp -^X T (G~ 1 + 2U)X j 



M 



det(I + 2GU) M / 2 



, (26) 



(27r)"/Vdet(G) 

where X = (X\, . . . ,X n ) T , and Uij = u^ij. This result is similar to that obtained for the 
Wishart matrix distribution [2T]. The next step is now to calculate the derivative of Z n with 
respect to «&. We set A = I + 2GU, and make use of the Jacobi formula 



d(det(A)) = det(A)Tr(A- 1 dA), 



(27) 



where dA is the variation of the matrix A when Uk undergoes a variation dut- From there 
one obtains 



dZ M {l-A- k l) 



Z n (u). 



(28) 



du k v ' 2 u k 

The second derivative can be calculated upon noting that (A'^kk = Ckk/ det(A), where Ckk 
is the k th element on the diagonal of the cofactor matrix of A. By construction, Ckk does 
not depend on Uk, and the following result is readily obtained 



d 2 Z n 
du\ 



u 



1 + 



M 



Mil- A~l) 



-i 2 



U k 



Z n (u). 



(29) 



As a check on these results, (uk) and (z/|) can be recomputed by setting u = respectively 
in Eq. ([281) and Eq. ([291. For small u, A" 1 



I 



2GU (first term of the Neumann series), 
so that (i/ fc ) = MG kk and (z/|) = 2MG 2 kk + (MG kk ) 2 are recovered. 
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C. Neutral density and ionization source 



The results obtained in the previous section, namely Eqs. (26 ) and (28 ), provide analytical 
expressions for the average neutral density and ionization source when the flux at the wall 
To is non-stochastic (slow recycling). The p th moment of the density is given by 

(N' M ) = (*)\W - (%)" det(I + 2 ^ GUt)M/ , . (30) 

where {Uk)ij = UkSij = e/v 5i<kSij- This result can be validated by investigating several 
limiting cases (see Appendix B). The average density (p = 1) is plotted on Fig. [2] as a 
function of the distance to the source in units of the mean free path I, for different values of 
the ratio a = A/7. As an example, for N e = 5 x 10 18 m -3 , T e = 50 eV and Eq = mv\j2 = 20 
eV (m is the mass of the atom), the ionization mean free path I is of the order of 20 cm for a 
D atom. In the same conditions, for a D2 molecule with Eq = 0.03 eV, I ~ 8 mm. Fig. |2]-a) 
shows the average density in the A — > 00 case, for different values of the fluctuation rate R. It 
is seen that the average density profiles decrease more slowly than the fluctuation free profile 



(R = 0, solid line), in accordance with Eq. (22), but that significant deviations occur only 



for large fluctuation rates above 50 %. Fig. 0-b) shows that as A tends to zero, the average 
density profile comes closer and closer to the turbulence free profile, as expected from the 
arguments developed in section A. The slowest decaying profile corresponds toa-> 00. The 
profiles for finite values of a are close to this limiting case provided iC A. On the size L = 8 I 
of the box, a = 90 provides an excellent approximation to the a = 00 case. The standard 



deviation of the density is obtained from Eq. (30) by (AN(xk)} = ((N (#&)) — (N(xk)) ) 



T2f„ \\ IAT(„ A\2\l/2 



and is plotted on Fig. [3]-a) for a = 90 and b) for a = 0.45 and M = 3, both as error bars on 
the average density profile and in the inset. The latter represents the typical dispersion of 
the different realizations of the density profile, and decreases with A in accordance with the 
results of section A (ergodic theorem). The standard deviation is zero for x = 0, because 
the boundary condition prescribed here is non stochastic. The average ionization source and 
its second moment are respectively given by 

r Bz 

(S(x k )) = -_ °__(U k ) = w(x k )(N(x k )), (31) 

v du k 
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FIG. 2: a) Average density profiles obtained with a non stochastic boundary condition, as a function 
of the distance of the source in units of mean free path I, for different values of the fluctuation rate. 
All the profiles are calculated for a = +oo. The solid line is the turbulence free profile, the dash 
dotted line corresponds to R=32%, the dashed line to R=58%, and the dotted line to R=82%. 
As R increases, the decay of the average density profiles becomes slower and slower, b) Average 
density profile obtained for R=82 %, for different values of the ratio of the integral scale to the 
mean free path a. The solid line is the turbulent free profile, the short-dotted line corresponds to 
a=0.45, the dashed line to a=2.7, the dash-dotted line to a=9 and the dotted line to a=+oo. The 
larger a, the slower the neutral particle density decays. On both plots, open circles represent the 
results of the Monte Carlo simulations carried out in section IIVI 



<S 2 (^)> 



v 



2 s<-> 



^)^(x k )(N"(x k )), 



(32) 



with 



w(x k ) 



Mv (1 - (A-, 1 ] 



kkl 



(33) 



2 e 

where we recall that A^ = I + 2GU&. vj(x k ) is an effective ionization rate, which depends 
on space. Its asymptotic behavior for large values of x k is obtained from Szego's theorem on 
Toeplitz matrices J2Z] in AppendixO In the continuum limit (e — > 0), and for the correlation 



function given by Eq. (14), we have 



Woo = lim w(x k ) 

k— >+oo 



[1 + 2(1/) 



T 



1/2 



(34) 
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a) 




b) 



FIG. 3: a) Plot of the average density profile as a function of the distance to the wall in units of 
mean free path 1, for R=82% and a = 90. The error bars are such that (N) ± (AN), where (AN) 
is the standard deviation profile, which represents the typical statistical dispersion of the density 
profiles. (AN) is plotted in the inset, normalized to the value of the density at the wall. For x = 0, 
(AN) = because of the non stochastic boundary condition, b) Same as a), but for a=0.45. The 
dispersion is clearly lower by a factor of about 2 compared to a). This result is consistent with 
what is expected from the ergodic theorem, which implies (AN) — > as a — > 0. On both plots, 
open circles represent the results of the Monte Carlo simulations carried out in section |IV| 



where the time r is such that r = 4A/Mt> . Eq. (34) shows that (S) and (N) have the same 



asymptotic decay, except in the infinite correlation length case where Wk, tends to 0. The 
effective ionization rate w is plotted on Fig. |4j for a = 0.45 and a = 900, together with the 
corresponding asymptotic values from Szego's theorem. This value is reached at a distance of 
a few correlation lengths from the wall, where w(xi) = (v). In the a = 0.45 case, e/A = 0.2, 
and the difference between the continuum limit and the discrete result given in appendix 



C (Eq. (C6)) is of the order of 2%. The average ionization source (S(xk)) is plotted on 
Fig. [5] for a = 0.45 and a = 900. The differences between the turbulence free and the 
averaged ionization source are much smaller than in the case of the neutral particle density, 
because the enhanced neutral particle penetration originating from low density realizations 
is compensated by the corresponding low ionization rates. The standard deviation is this 
time maximum at the wall, since the fluctuations of S(0) directly reflects those of iV e (0). 



This can be checked directly by forming (S (0)) — (S'(O)) from Eq. (31) and (32). Finally, 
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FIG. 4: Plot of the effective ionization rate w in units of (u), as a function of the distance to the 
source for two values of the integral scale, namely a=0.45 and a=900 (in units of 1). For a=0.45, 
w(x) decays towards a constant value of the order of 0.8, in accordance with Szego's theorem 



prediction, Eq. (C6). In practice, this value is reached after a distance of the order of twice the 
integral scale. 

we come back to the neutral particle density, this time when the boundary condition T is 
stochastic (fast recycling, denoted by the F superscript). In our model, the recycling flux 
is directly proportional to the plasma density, i.e. to the ionization rate, and the average 
neutral particle density and its second moment are given by 



(N F (x k )) 



(N F (x k ) 2 ) 



(To) dZ 



v Q {v) du } 



(Ui 



™{Xk) 



(N(x k )) 



M\ / w(x k ) 



1 + T 



(*> 



(N 2 (x k )) 



(35) 



(36) 



according to Eq. (24), (28), and (29) and using (A^ 1 )^ = (A^ 1 )^ (see Appendix O. 
Therefore, the results are the same as for the ionization source, up to a dimensional factor. 
In this case, the source strength plays the same role as the ionization rate for S in the slow 
recycling case. 



The results presented in this section highlight the essential role played by the ratio a of 
the turbulence correlation length to the neutral mean free path. This provides an important 
clue to understand 2D and scattering effects. Moreover, if the monokinetic approximation 
is relaxed, linearity implies that the average density can be obtained by further averaging 
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FIG. 5: a) Plot of the average ionization source (S) as a function of the distance to the wall in 
units of the mean free path, for R = 82%. The dashed line is the turbulence free case, the solid line 
corresponds to a = 0.45, while for the dotted line a = 90. The differences between the averaged 
and the turbulent free profiles are smaller than for the average density, b) Plot of the average 
ionization source, for R=82% and a=1.3, with errors bars representing the standard deviation 
(AS). In contrast with the neutral particle density, the dispersion is the largest at the wall. Open 



circles represent the results of the Monte Carlo simulations carried out in section IV 



Eq. (30) over the initial velocity t> distribution function. This is equivalent to average over 



a, since the neutral mean free path is directly proportional to Vq, so that in the end one 
might have properly scaled contributions from both a <C 1 and a>l regimes. For instance, 
atoms created by molecular dissociation in edge plasmas can have energies as low as 0.2 eV, 
while backscattering at the wall leads to atoms having energies of several hundreds of eV. 

D. Influence of the correlation function shape 

In the previous section, the role of the integral scale A has been discussed for a correlation 



function C(r) of the form given by Eq. (14). However, Eq. (30) is valid for any positive 
(a restriction of the multivariate Gamma model) C(r). In this section, we address the 
sensitivity of the average density profile to the actual expression of C(r), for the same 
value of the integral scale A = J C(r)dr/C(0). We first consider the following family of 
correlation functions C q (r) 
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C q (r) 



a 



(l + (r/r ) 2 )«' (37) 

for which r = (2A/7r)(2g — 2)!!/(2g — 3)!!. In addition, we consider a Gaussian correlation 



function C g (r), which decays faster than Eq. (14). These various functions are compared 
on Fig. [6] a), with q = 1 and q = 10 for C q , for the same value of a = 1.3. The resulting 
density profiles are plotted on Fig. [6]b), and are found very close from each other. This 
observation is found to be valid for any value of a, the largest deviations being found for 
a — ¥ +oo, where the effect of fluctuations is the largest (for a — > all the density profiles 
tend to the fluctuation free result). Therefore, for a given integral scale, the shape of the 
correlation function is likely to have only a minor effect on the average quantities. In other 
words, the main control parameter is the integral scale of the fluctuations. 
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FIG. 6: a) Plot of an exponential (Eq. (14), solid line), gaussian (dotted line), and power law (Eq. 



(37) for q=2, dashed line, and for q=10, dash-dotted line) correlation functions, with the same 



integral scale (i.e. the same area under the curve) such that a=1.3. b) Plot of the corresponding 
average density profiles as a function of the distance to the source in units of the mean free path. 
The different profiles are found to be very similar. 



E. 2D effects 

We now investigate 2D effects, since Monte Carlo calculations including scattering (i.e. 
charge exchange) will be performed in 2D. The neutral density profile is expected to be more 
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peaked towards the wall in 2D than in ID, because neutrals which have large v y velocity 
components will be ionized closer from the wall (along x) than those having small or zero 
v y component. In fact, if 9 is the angle between the initial velocity {i.e. f2) and (Ox), 
x = s/ cos 8. Assuming a cosine distribution for the particle source at the wall, the following 
solution for / in an homogeneous plasma is readily obtained by averaging Eq. (|6| in the 
scattering free limit over angles 

f(x,v) = 2N 8(v -v )—E 2 (-) , (38) 

where E2 is an exponential integral [26j. The following result for the average density is easily 



obtained from Eq. (30) in the slow recycling case 



where \x = cos 9. The latter expression of (N 2 D{xk)) would be strictly exact in the continuum 
limit, since in principle G and U^ depend on /i even if isotropy (namely C(r— r') = C(|r— r'|)) 



is assumed, because of finite grid step effects. However, Eq. (39) provides very good results 
in practice for e/C = 10~ 2 . The density profiles calculated for M = 3 and a = 00 are 
compared to their ID counterparts on Fig. [7} Both 2D density profiles are found to be 
more peaked towards the wall than in ID. It should be noted that the angle average can be 
interpreted as an initial velocity average, with a velocity distribution fo(v) such that / = 
for v > Vq. 

IV. IMPLEMENTATION IN THE EIRENE MONTE CARLO CODE 

We now discuss the numerical implementation of our stochastic model in the parallelized 
(MPI) EIRENE Monte Carlo solver [5]. We have adopted the following procedure. The 
plasma parameters, say the electron density, are generated by sampling n 2 values from 



the multivariate Gamma distribution defined in section II C , then Boltzmann's equation is 
solved in this plasma background. These steps are repeated until statistical noise due to 
random sampling of the background medium is reduced to a negligible level, compared to the 
intrinsic Monte Carlo noise for each realization of the stochastic plasma fields. The numerical 
implementation is validated against the analytical results obtained in the previous section. 
The effects of scattering (i.e. charge exchange in our context) are then finally investigated. 
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FIG. 7: Comparison between the 2D and the ID cases for M = 3 and a = oo. The solid line 
and the dash-dotted lines are the turbulence free profiles in ID and 2D respectively. The dotted 
and the dashed-double dotted line are the average density profiles in ID and 2D respectively. The 
effect of fluctuations is similar in both cases. Open circles represent the results of the Monte Carlo 
simulations carried out in section ITVT . 



A. Sampling from the multivariate Gamma distribution 



From the definition of the multivariate Gamma distribution given in II C it is clear that 
sampling from the latter requires generating M series of n Gaussian numbers with correlation 
matrix G. By construction, G is symmetric and positive definite so that it can be Cholesky 
factorized in terms of a lower triangular matrix L such that G = LL T (where superscript 
T denotes matrix transpose) [28] ■ Once L is calculated, and if Yx, . . . , Y n are n independent 
Gaussian numbers, then X = LY is a vector of n elements Xi which have correlation matrix 
G. This stems from the fact that 



X T G- X X = (L- 1 X) T L- 1 X = Y T Y, (40) 

so that by the law of transformation of probability densities, if Y has a multivariate normal 



law, X = LY is distributed according to Eq. (11). We present on Fig. 8 a 2D map 



(100 x 100 cells) of variables Ui, obtained from Eq. (12) using the same sample of Y. The 



matrix G (size 10 4 x 10 4 ) is constructed from Eq. (15), with M = 3 (fluctuation rate 
R ~ 82 %) and for values of X/C from 10~ 3 to 2. The same color scale is used for all four 
sub-figures. The same underlying pattern can be recognized, but is progressively smeared 
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out as A increases. For values of A large compared to the size of the domain, the field v{r) 
becomes constant in space (but the value taken by v is different for each realization). In 
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d) 

FIG. 8: 100 x 100 2D map sampled from the multivariate Gamma distribution, for M = 3 and a) 
X/C = 10~ 3 , b) X/C = 2.5 x 10~ 2 , c) X/C = 0.15 d) X/C = 2. Note that the average density is 
(JV e ) = 5 x 10 12 cm" 3 . 



practice, straightforward application of this technique in 2D is currently limited to a domain 
size of the order M = 10 5 cells, mainly because of memory issues. In fact, for n = 300, the 
correlation matrix has 4 x 10 9 independent elements, i.e. ~ 30 Gb in single precision. 
Cholesky factorization requires A/" 3 /3 operations for a jV x J\f matrix, that is about 100 s 
for M = 300 x 300 on a 80 GFlop workstation. As a result, we are effectively restricted in 
2D to X/C > 3 x 10 -2 . This is clear on Fig. p^J where density is constant on a cell of size 
£/100, whatever the value of A. As a result, choosing X/C = 10 -3 in this case ensures that 
two neighboring cells are uncorrelated, but the effective correlation length is of order £/100. 
In ID, if there are n cells, the correlation matrix isnxn and X/C > 10~ 6 . 
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B. Source of neutrals 



We discuss here the implementation of the stochastic boundary condition introduced 



in section II D where the magnitude of the neutral source is proportional to the plasma 



density at the wall (fast recycling). When the plasma density background has been sampled, 
r p (y) = N e (0,y)Vuob is computed and then used, with the proper normalization, as a PDF 
for creating neutrals along the y direction, at position x — 0. The results are finally rescaled 
using for each realization the total incoming flux by integrating T p over y This procedure 
is illustrated by plotting a particular realization of the density map on Fig. M a), and the 
corresponding PDF on Fig. Mb). 
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FIG. 9: a) sampled plasma density map for A/£ = 0.15, the wall is at x = b) corresponding 
neutral birth pdf, plotted as a function of cell number (cell number 100 corresponds to y = C). 



C. Validation 

The implementation of our stochastic model for the plasma background in the EIRENE 
Monte Carlo solver has been checked by running cases corresponding to the model described 
in the previous sections. Charge exchange and the density dependence of the rate coefficient 
o~i v e are turned off. To simulate a ID (scattering free) problem, a point source is used in a 
2D slab, and all the neutrals are launched with the same velocity v$ along the x axis. In this 
extremely simplified case, the conditional expectation estimator implemented in EIRENE 
provides a zero variance Monte Carlo scheme, tracking only one particle |29|. Numerical 
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results for the average density and the standard deviation after many realizations of the 
stochastic background are plotted on Fig. [2] (open circles) for a = oo with M = 3, 6 and 
20, and for M = 3 with a = 0.45, 2.7, 9 and oo. In the simulations, neutrals are removed 
(absorbed) when crossing the x = £ surface. The standard deviations are also identical 
to those obtained analytically, as shown on Fig. [3j again with open circles. The same 
agreement is obtained for the average ionization source, and the corresponding standard 



deviation (see Fig. O. The fast recycling case has been validated against Eq. (35) and 



(36). Finally, numerical results in 2D are validated against Eq. (39) on Fig. [7| For these 
runs, periodic boundary conditions are imposed at y = and y — C. This exercise provides 
strong cross checks between the numerical and the analytical results, and allows estimating 
the number of realizations of the plasma background needed to obtain reliable results for 
the averages and the standard deviations of the various quantities of interest. In fact, in 
the worst cases, i.e. for large values of A, 10 4 realizations provide very good estimates for 
the two first moments, but these are still not enough for higher order moments (which can 



be calculated for the density in the slow recycling case from Eq. (30)). In practice, for 
small values of A, less iterations are needed because the density profiles along x are space 
averages along the y direction, and that different stripes defined by yi < y < yi + \ such that 
yi + i — yi > A can in a first approximation be seen as independent realizations. 



D. Calculations including scattering 

In this last section, we investigate whether taking scattering into account changes the 
overall picture that was obtained in the 2D scattering free case. We limit ourselves to the 
effect of density fluctuations. In the case of D atoms in a plasma, this amounts to neglect 
both ion temperature and velocity fluctuations, since the properties of charge exchange 
(CX) neutrals depend on those of the background ion velocity distribution. In the edge 
plasmas of tokamaks, a simple estimate of the turbulent velocities shows that the latter are 
small compared to the thermal velocity [30] , so that the velocity of CX neutrals should not 
be strongly affected by turbulence. In an homogeneous background medium, if scattering 
conserves kinetic energy (one speed transport model), the neutral decay length should be 
smaller than in the scattering free case because trajectories are no longer ballistic. In the 
limit of vanishing scattering mean free path (compared to plasma inhomogeneities), the 
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transport becomes diffusive. If the initial neutral velocity t>o differs from the thermal ion 
velocity v t h, neutrals which have undergone charge exchange have different mean free paths 
than first generation neutrals. The latter have a parameter a calculated from E = l/2mt>Q, 
while all subsequent generations have a value a calculated from 3T,/2. To summarize, when 
scattering is retained, at least two more parameters enter the problem. The first one is 
the ratio of the scattering to the absorption rate, b = a s v/a a v, where a a and a s are the 
corresponding cross sections. The second parameter, denoted by d, is the ratio of the energy 
of scattered neutrals to their initial energy, namely d = 3Ti/2E (in the one speed problem, 
the velocity after charge exchange is calculated from v cx = ^3kTi/m so as to ensure proper 
thermalization). It should be noted that b and d are closely related to the parameters which 
control the analytical solution of Boltzmann's equation obtained in simplified ID geometry 
with mono-kinetic ions by Smirnov, namely /3 = a s v/(a a v + a s v) and uq, the neutral initial 
velocity in units of ion thermal velocity [31]. In the calculations presented below, periodic 
boundary conditions are implemented for y = and y = C, the x = surface is reflecting 



and x = C absorbing. The effect of b is investigated on figures 10 a) and b) for d = 1, where 
results respectively obtained for b = 0.75 and b = 3.5 are compared to those obtained for 
the scattering free case (6 = 0). The calculations confirm that the neutral density is larger 
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FIG. 10: a) Density profiles calculated in the slow recycling case for A = oo, b = 0.75, d = 1, in 
the fluctuation free case (solid line) and for R = 82% (dotted line), compared to the scattering free 
results (b = 0, fluctuation free (dotted line) and R = 82% (dash-dotted line)), b) Same as a), but 
for b = 3.5. 
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FIG. 11: a) Density profiles calculated in the slow recycling case for A = oo, 6 = 3. 5, d = 0.1, in 
the fluctuation free case (solid line) and for R = 82% (dotted line), compared to the d = 1 case, 
fluctuation free (dotted line) and R = 82% (dash-dotted line), b) Same as a), but for d = 10. 



close to the wall when scattering is included, and show that the effects of fluctuations are 
similar than in the scattering free case, even for b = 3.5. The effect of the parameter d is 



shown on Fig. 11, where density profiles are plotted for b = 3.5, d = 0.1 and d = 10, and 
compared to the d = 1 case (for which the neutral particle energy does not change in the 
scattering event). In the latter case, the effect of fluctuations is clearly reduced compared 
to d — 1 or d — 0.1, as expected. In fact, after charge exchange, the neutral particle have 
a higher energy, hence the ratio a of the turbulent integral scale to the mean free path is 
smaller. 



V. CONCLUSIONS AND PERSPECTIVES 

We have presented a stochastic model which allows investigating linear transport in a 
turbulent background. The statistical properties of turbulence have been described by a 
multivariate Gamma law, because the latter has three major attractive features. First, aver- 
ages can be calculated analytically in the scattering free case, which is relevant for molecules 
and sputtered neutral impurities in tokamaks. Next, it is conveniently implemented in a 
Monte Carlo code, here the EIRENE Monte Carlo Solver, so that analytical and numerical 
results can be cross-checked. Finally, the univariate Gamma distribution [i.e. the one 
point marginal) provides a good description of plasma density fluctuations in the outer 
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edge of tokamaks. The multivariate Gamma law is fully specified by its correlation matrix, 
which is closely related to the spatial spectra of fluctuations. In the frame of this stochastic 
model, the average neutral density and ionization source profiles have been investigated. In 
the scattering free case, two main control parameters have been identified. The first one 
is related to the properties of the source, namely whether it is stochastic or not. In the 
frame of neutral transport in plasmas, this is controlled by the ordering between recycling 
and turbulence time scales. We thus distinguish two limiting cases, namely slow and fast 
recycling, the relevance of which depends both on the elementary process at play and the 
wall status (i.e. saturated in hydrogen or not). The second control parameter is the ratio 
a of the turbulence correlation length (or integral scale) to the average neutral mean free 
path. For small values of a, the ergodic theorem implies that the average quantities can be 
calculated by replacing the absorption rate by its average in Boltzmann's equation, thus 
recovering the so-called atomic- mix limit in Ref. [TJ. This does not mean that turbulence 
has no effect, since the average ionization rate can differ from the ionization rate calculated 
for the average plasma parameters, in particular when temperature fluctuations are present. 
For non-stochastic boundary conditions, fluctuations enhance neutral penetration, in the 
sense that the decay of the average density profile is slower than in the turbulence free case. 
This effect becomes more and more pronounced as a is increased. The differences between 
the average ionization source, and the ionization source calculated in the turbulent free case 
are less significant. The physical reason for this is clear. Realizations which have typically 
low density lead to deep neutral penetration, but at the same time to low values for the 
ionization rate. In other words, the source, being homogeneous in space, is significant even 
in regions where the plasma has low density. In contrast, in the fast recycling case the 
source is concentrated where the density is high, so that neutrals are more likely to be 
ionized close to the wall provided the electron temperature is high enough. In this case, the 
ionization source becomes more peaked towards the wall. These conclusions remain valid 
in the presence of scattering, as shown by numerical calculations. They provide a clear 
physical picture of the role of fluctuations on neutral transport, and are therefore expected 
to remain qualitatively correct when fluctuations obey a different statistics. The effect of 
turbulence is significant only for large fluctuation rates (> 50 %). In the fusion context, 
this means that our results are of interest mainly for SOL plasmas. The confined plasma 
might nevertheless be affected, through a reduction of the SOL screening efficiency. 
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The work presented here is a first step towards evaluating the importance of turbulence on 
neutral particle transport in fusion plasmas. Further work will now focus on implementing 
a more realistic description of the outer edge of the plasma, i. e. to relax both the statistical 
homogeneity and stationarity assumption. The role of temperature fluctuations must also 
be thoroughly studied, because of threshold effects in the ionization cross section. Both can 
be addressed in the frame of our stochastic model, or rely on the output of a turbulence 
code. A direct coupling to such a code would finally allow to relax the passive assumption, 
i.e. to investigate the back reaction of neutrals on turbulence. 
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Appendix A: Expression of the chi-2 PDF 



We consider here the one point marginal W\ of Wn defined by Eq. (17). The integration 
of Wn over n — 1 variables is straightforward, and leads to 



m 




W 1 (u)= dX x ... dX M exp -^ i h~E X n- ( A1 ) 



2a 2 
To calculate W(v) we note that the integral is of the form 



W x {u) = I dX 1 ... I dX M P(X u . . . , X M )5 (fiX,, . . . , X M )) . (A2) 

Using the co-area formula, we get 
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W 1 {y)= [ dS(X) T - i —P(X 1 ,...,X M ), (A3) 

Ju I V J\ 

where U is the set of point in R M such that /(X) = 0, and S a measure on this set. Now, 
all the terms in the integral involve only J^ i=1 A 2 an d can simply be taken out, hence 

Wi{y) = —pL ^exp ( — ^-) [ dS(X). (A4) 

The remaining integral is the area of the hypersphere of radius r = z/ 1//2 in dimension M, so 

that 

^m- ^^m) ^" «»(-£)■ (A5) 

which is the usual expression of the chi-2 distribution with M degrees of freedom and scale 
parameter a ([32], P- 452). The later is nothing else but a Gamma distribution with a shape 
factor a = M/2, a scale factor 2cr 2 and a displacement 7 = ([32], P- 337). 

Appendix B: Limiting cases for the average density profiles 



We investigate here three limiting cases of Eq. (30) with p — 1 and of Eq. (35), namely 
the small fluctuation limit, the infinite correlation length limit (A — > oo), and its opposite 
A — y 0. First, the vanishing fluctuations limit corresponds to M — > +oo and Gij oc \jM — > 
(i.e. a = z^oa/2/M), so that (u) —¥ Vq (where Vq is the ionization rate in the turbulence free 
case) and ((^ 2 )) — > 0. As M becomes large, 

det(I + 2GU fe ) ~ 1 + 2Tr(GU fc ) = 1 + ^, (Bl) 



so that using Xk = ke Eq. (30) reduces to 



(N(x k )} ~ — exp-is x k /v , (B2) 



as it should. In the case of Eq. (35), the limit of w must also be studied. For small 
fluctuations, we have A -1 ~ I — 2GUfc, hence 1 — (^4 _1 )« — 2e(u)/MvoSi<k, so that w 
reduces to (i>). Now we consider the A — > limit, where G tends to a diagonal matrix 
Grf such that (Gd)u = a/y2M. The calculation of the determinant is then trivial, and we 
obtain 
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(N (x k )) 



1 



1 M/2 



(B3) 



v l(l + 2x k G /kv ) k _ 

where G = (v)/M. But when A — > 0, we must also impose e — > 0, k — > oo, fee = x, so that 
the average density profile reduces to 



(N (x)) ~ — exp-(i/)a;/-Uo, 

wo 



(B4) 



as expected from the ergodic theorem. In the the fast recycling case, it is straightforward to 



show that w — > (u), so that Eq. (35) also reduces to Eq. (B4). Finally, we study the infinite 
correlation length limit A — > oo. In this limit, GU& becomes proportional to an x n matrix 
Bfc such that its elements are equal to (-Bfe)y = 8j< k - This matrix has two eigenvalues, 
namely with rank n-1 and k with rank 1. The calculation of the determinant is then again 
trivial, and leads to 



(N^Xk)) 



Ell 



M/2 



(B5) 



l + 2x k G /v _ 

For large correlation lengths, the plasma background becomes uniform, so that this result 
can be recovered from a simple integral, namely 



(Woo(aO) 



^o Jo 



+oo 



dvW\{y} exp(—vx/vo), 



(B6) 



where W\ is the univariate chi-2 distribution with M degrees of freedom given by Eq. ( Al ) 



Furthermore, for A — > +oo, A k can be expressed as A k = I + (2Goe/vo)Bfc. It is easily 
shown by induction that (B^)™ = k n ~ 1 'Q k) so that expanding A^ 1 in Neumann series leads 

to 



n=0 V / 



/ - 2 J^1 (l + 2G ° Xk 



-1 



B, 



(B7) 



vq V v o 

where x k -- ek. In the continuum limit where e — > 0, k — > +oo and x k = x, the effective 
ionization rate w reduces to 



w = (v) 1 + 



2G x k \ x 

Vq ) 



(B8) 



This result is consistent with that obtained from the integral 
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T-l ^ + 00 

(N^ix)) = — °- / dvW 1 (u)vexp(-vx/v ). (B9) 

v o\ u ) Jo 

Appendix C: Asymptotic behavior of w(x), an application of Szego's theorem 

Consider a set of symmetric Toeplitz matrices T n (n x n) with n = l,...,+oo, i.e. 
matrices such that (T n )ij = (T R )i_ij-_i for i > 2. A Toeplitz matrix is completely specified 
by 2n — 1 entries t m , with m G \—n — 1, n + 1], where positive values of m refer to the upper 
part of the matrix. For a symmetric matrix t_ m = t m . We define a function f(p) by 

+00 

m=— 00 

with pe [0, 27r]. Szego's theorem implies [27] 

llm !?!%I = * f*ln/ W = L„. (02) 



,^-+00 det T n 27r jq 

We apply this result to a set of k x fc Toeplitz matrices T& such that T^ = !& + lejv^G^. In 
this case 



T/26 

t m = <5 m o + C exp - — , (C3) 



where Go = cr/y2M. It is easily checked that det(Tfc) = det(Ak) (where we recall that 
Ak = I n + 2GUk is a n x n matrix ). The matrix elements (A^ 1 )^ can be written in terms 
of the corresponding cof actors Cu, as 

(A- k % = (-If -^ = d ^^ 1 , (C4) 

det(T fc ) det(T fc ) 

where the last equality holds for i = 1 and z = k only. Therefore, for large k, (A' k 1 )u ~ L^. 

Furthermore, (A^ 1 )].! = (A^T 1 )*^ because the corresponding minors, obtained by removing 

either the first or the last line and column of the Ak matrix, are equal. A straightforward 



application of Eq. (CI) with the t m defined by Eq. (C3) leads to 



f(p) = 1 + -G tanh (±) ) (C5) 

v V2A/ 1 — cos(p)/ cosh(e/2A) 

The integral in Eq. (tC2J) can be calculated analytically, yielding 
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;i + tanh(e/2A)) 



g + (g*-l/cOSh 2 (e/2\)) 



1/2' 



where 



(C6) 



9 = H G tanh ( — J . 

fn V2A/ 



Finally, in the continuous limit where e — > (in practice e/A <C 1) 



,. 1 — -^oo 

hm 



1 

2A 



from which Eq. (34) is obtained 



1 + 



oGrnA 



v 



1/2 



(C7) 



(C8) 
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